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ABSTRACT 

We construct equilibrium configurations of magnetized, two-fluid neutron stars us¬ 
ing an iterative numerical method. Working in Newtonian framework we assume that 
the neutron star has two regions: the core, which is modelled as a two-component 
fluid consisting of type-II superconducting protons and superfluid neutrons, and the 
crust, a region composed of normal matter. Taking a new step towards more complete 
equilibrium models, we include the effect of entrainment, which implies that a mag¬ 
netic force acts on neutrons, too. We consider purely poloidal field cases and present 
improvements to an earlier numerical scheme for solving equilibrium equations, by 
introducing new convergence criteria. We find that entrainment results in qualitative 
differences in the structure of field lines along the magnetic axis. 
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1 INTRODUCTION 


Neutron stars are ideal astrophysical laboratories. Their 
extreme properties demand an interdisciplinary approach, 
while most of the existing models attempt to synthesize the 
various features that neutron stars are predicted to possess. 
Of particular interest is the construction of numerical so¬ 
lutions for equilibrium configurations of magnetized neu¬ 
tron stars with realistic assumptions about their interior 
structure. The Hachisu self-consistent field method (HSCF) 
(H achisu] l986) is the starting point of many modern meth¬ 
ods for constructing (initially non-magnetized) equilibria. 
This method was extended by Tomimura & Eriguchi (2005) 


so that magnetic fields with a strong poloidal and a weak 


toroidal component were included. Lander et al. (20121; Lan 


der (2013) and Lander (20141 further extended the methods, 


using a two-fluid description, including type-II superconduc¬ 
tivity for the neutron star core. In these latter models, the 
star is divided in two regions, the superconducting core and 
the normal infinitely conducting crust. 

Here, we also take into account, for the first time in such 
numerical models, the effect of entrainment. This comes by 
assuming an interaction between the momenta of the two 
fluids, resulting in an additional force acting on neutrons, 
see IGlampcdak fs et al.|(2011[ ) and references therein. 

Furthermore, because the number of equations to be 
solved has increased, we show that the simple convergence 
criteria introduced in|HachIsu|( 1986| are no longer adequate. 
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We present new convergence criteria, related to fundamental 
quantities of the model and show how one obtains more 
accurate numerical solutions. 

Although our code is built in a way that it is possi¬ 
ble to numerically evaluate both purely poloidal and mixed 
poloidal-toroidal equilibrium models, we focus here on the 
case of purely poloidal fields. The reason for doing so is that 
in these models B$ <g B po i. We compute models with en¬ 
trainment for a magnetic field strength of ~ 10 15 G at the 
pole and compare to corresponding models without entrain¬ 
ment, that were previously obtained in Lander (20141. 

In the following sections we first briefly discuss entrain¬ 
ment, then we describe the model. We demonstrate our new 
convergence criteria along with the numerical method and 
finally we provide results for equilibria with and without 
entrainment, and compare these. 


2 ENTRAINMENT IN NEUTRON STARS 

It is theorized that neutron stars exhibit entrainment, a 
phenomenon occurring in the case of mixed fluids. Due to 
entrainment, the momenta of the components are coupled, 
which means that the momentum of each component is de¬ 
pendent on both fluid velocities. Thus, if a component is 
flowing, the other component will flow as well. Following 
|Prix et al.| ( |2002[ ) and |Glampedakis et al.| ( |2011| ) the en¬ 
trainment can be simply described through a dimensionless 
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function e* defined by 


£* 


1 £n £p 
1 £r\ 


(i) 


where e n and e p are entrainment coefficients, satisfying, by 
definition, the relation p p e p = p n £ n , where p p and p n are 
the proton and neutron mass densities. We here consider 
the slow-rotation approximation, neglecting rotation related 
terms in the equations of hydrostationary equilibrium as the 
shape of the star is still spherical.The coefficient £ p can be 
written in terms of the bare proton mass m p and the ef¬ 
fective mass m* acquired by each proton as a result of the 
entrainment, as 



and thus £* can be written as 

_ p p (6,, — 1) + PiA 
p p (5* - 1) + p n 


( 2 ) 

( 3 ) 


where we have set <5* := yyf -, which depends on the den¬ 
sity and thus on the equation of state. In general, <5* varies 
depending on the model taken into account. Some of these 
models are discussed in |ChameT| p008). Here, we assume 
that it is constant throughout the core. Notice that in the 
case that m p = m p , e* = 1 and the entrainment effect van¬ 
ishes. 

Entrainment allows neutrons to experience a magnetic 
force. Normally, this force would be absent since neutrons 
are neutral and thus do not interact with the magnetic field. 
However, since the protons are assumed to be superconduct¬ 
ing and the neutrons superfluid, the drag between the proton 
and neutron fluid components provides a mechanism of in¬ 
teraction between neutrons and the magnetic field, which is 


enacted through the entrainment force on neutrons (Alpar 
et al. |1984 1 . 


3 THE MODEL 

3.1 General description 

Working in Newtonian framework we assume a magnetized, 
axisymmetric neutron star in a stationary state. The axis of 
symmetry is parallel to the magnetic dipole axis. The star 
is assumed to have only two regions: the core and the crust. 
The core, a region that starts at the center and extends to 
roughly 90% of the radius of the star (the crust-core bound¬ 
ary) , is assumed to consist of superfluid neutrons and type-II 
superconducting protons, where the protons are subject to a 
magnetic force different from the normal ideal MHD caseQ 
while neutrons will also be subject to a magnetic force due 
to entrainment. The crust, which lies between the crust-core 
boundary and the surface of the star, is composed of normal 
matter, which is subject to the standard Lorentz force. The 
crust is assumed to be in a relaxed state, without tangential 
stresses and can thus be described as a single fluid. Outside 


1 For simplicity, we ignore the presence of electrons, as well as 
more exotic particle species that may appear at high densities in 
the core. 


the star we assume vacuum. For purely numerical conve¬ 
nience we adopt the crust model described in Lander (20131, 


Lander (20141 and Lander et al. (20121. In the future, using 


realistic equations of state, it should be possible to adopt 
more realistic proton and neutron fractions throughout the 
star. However, for evaluating the entrainment effects in the 
core, the assumptions made for the crust are not crucial. 

We express the main equations in cylindrical coordi¬ 
nates ( U7,(p,z ), while in the numerical implementation we 
use spherical polar coordinates {r,9). Vector components 
are expressed with respect to the unit vectors (e^e^e^) 
while we provide some results in both dimensionless form 
and Gauss-CGS units. 


3.2 Two-fluid description 

In the core, we need to consider two separate equations for 
the hydrostationary equilibrium (obtained by setting the 
neutron and proton velocities equal to zero in the general 
MHD equations presented in Glampedakis et al.|2011 ): 

Fn 

Pn ’ 


V/i n + V<E> 

and 

VA P -(- V 4? 


( 4 ) 


Pp 


( 5 ) 


where An, Ap are the neutron and proton chemical potentials 


and 4? g is the gravitational potential (see Mendell 1991). 


F mag is the magnetic force acting on protons and F n is the 
magnetic force on neutrons due to entrainment. We will dis¬ 
cuss chemical potentials, as well as the form of the magnetic 
force in the following sections. 

The gravitational held is related to the total density p 
through 

V 2< hg = 4nGp = 47rG(p n + p p ). (6) 

We subtract 0 from 0 to find 
Fmag _ F 
Pp P. 

which we use along with |5| (instead of using S and Q). 
Since all quantities on the left side of Q and |5^ are gradi¬ 
ents of scalars, the same must hold for the right side, so 

F 


V(Ap - An) = 


( 7 ) 


= VM, 

Pp 

and 

F 

— = VJV, 

Pn 


( 8 ) 


( 9 ) 


where M, N are scalar functions. This result is useful, since 
( 5|7 1 then admit first integrals 

Ap T 4? g — Al — Cp = 0, (19) 

and 


Ap — An — Af -|- N — Cdif — 0, 


( 11 ) 


where C p and Cdif are some integration constants. The 
above first integrals will be used in the numerical implemen¬ 
tation. The magnetic field B needs to satisfy the divergence 
free constraint 


V ■ B = 0, 


( 12 ) 
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which implies that in axisymmetry it can be decomposed as 

B = —Vu x + B^ep, (13) 

zu 

where u is the streamfunction describing the magnetic field. 
The system is closed with an equation of state in terms of 
an energy functional, which will be discussed in the next 
section. 


3.3 Equation of state 

In the two-fluid description, equations @> 0 are formu¬ 
lated in terms of the chemical potentials p p , p n instead of 
the pressure P. We assume an equation of state in terms of 
the energy functional £(p n ,p p ) given by 

£ = k n pl + ^ + fc p p p +w? , (14) 

where k n ,k p are constants, and N n ,N p are the polytropic 
indices for neutrons and protons respectively. This is essen¬ 
tially an extension of the usual polytropic equation of state. 
The chemical potentials are defined through 

= fc n (l + ^)#, (15a) 

= k p (l+^jpf*. (15b) 

In the limit of a single-fluid description, one recovers the 
usual equations in terms of pressure, by 

VP = PnV/Xn + PpV/ip. (16) 



d£ 

Pp — 

dp P 

and 



d£ 

/^n = 

dp n 


3.4 Type-II superconducting core 


We assume type-II superconductivity for the protons in the 
core. Thus, the magnetic force is no longer the familiar 
Lorentz force, but it is replaced by a flux tube tension force 
(Glampcdakis et al.|2011 ) 


F 


mag 


1 

4-7T 


B x (V x Hoi) + PpV 



(17) 


where H c i = H c i(p p ,p n ) is the first critical field given by 
Hoi = H c iH and B is the unit tangent vector to the mag¬ 
netic field (B = B /B = B^e^ + B^e^ + B z e z ), with B 
the norm of the magnetic field. The norm of the first critical 
field is approximated by 


Hd{p P , Pn) = h c — , (18) 

where h c is some constant and £* is the entrainment function 
0- The neutron force F n is given by 

F " * -£ v • < I9 > 

This force on neutrons exists due to the coupling of proton 
and neutron flows at the vicinity of each fluxtube. These 
flows exist on a mesoscopic level and are different from the 
macroscopic fluid velocities which enter the fluid motion 


equations (see Glampedakis et al.||2011). 

We define D p and 

D n as 


D - dH * 

P ■" d Pp ’ 

(20a) 

and 


D - dH 

11 ' dp n ' 

(20b) 


If e* is constant, then H c i = H c i(p p ) and hence D n = 0. In 
that case, the force acting on neutrons vanishes. Combining 
|9]| and (191 it is obvious that 

N = ( 21 ) 

We extend the derivation for the type-II superconducting 
equivalent of the Grad-Shafranov equation for the magnetic 
field (L ander|2014 ) in the case that H c i = 7f c i(p p ,p n ) 

a . _ Vn-Vu _ 2 „ yjdy n2j d/ sc , 0 ^ 

A * u ~ n ppU du n/sc d«’ (22) 

where 

A* U ee _ 1 = v 2 (, (23) 

(24) 

(25) 


/ d 2 u 1 du d 2 u\ 

i - V 2 1 

f u sin</\ 

\ dzu 2 zu dzu ' dz 2 ) 

sin <j> 

l ™ ) 


while functions / sc and y(u) are defined through 
y(u) := 4nM sc + BD P , 


B, 


/sc(u) := ru^H cl . 

We note here that the contribution of the neutron force F n 
to the Grad-Shafranov equation is implicit. B is related to 
u through 

B = VB B = H c * - 


Cl V»c 2 l- 2 -/lc 

where we have set 

B Viil 


n := 


Ha \Jro 2 Hl i - / a 2 c 


(26) 


(27) 


We have denoted M and / with subscript sc in order to dis¬ 
tinguish them from their normal matter counterparts. The 
derivation of the Grad-Shafranov equation (221 is shown in 
section EH 


3.5 Normal matter crust 

In this section we discuss the form for the B field in the 
crust region. There, we assume normal, perfectly conducting 
matter and hence, the governing equations are those of the 
ideal MHD used in the first part. The magnetic force is 

F mag = (V x B) x B, (28) 

471 

while the neutron force vanishes 

F n = 0. (29) 

The Grad-Shafranov equation for this case is 

2 d M n . d/N , , 

A*u = -47TO7p p — - /n-i— , (30) 

dn dw 

where we have denoted the functions M, f with subscript 
N to distinguish them from their superconducting counter¬ 
parts. It can be shown that both Mn and /n are functions 
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of u. As in the superconducting case, /n is related with the 
toroidal part of the magnetic field B$ through 

/n(m) = vjBj,. (31) 


3.6 The exterior 

The exterior of the star is assumed to be a perfect vacuum, 
since we have not assumed the presence of a magnetosphere. 
Therefore, the matter density vanishes in the outer region 
and the equation governing the magnetic field is 

A*w = 0. (32) 

It is possible to model a magnetosphere by assuming that 
part of the toroidal component of the magnetic field exceeds 
the surface of the star, as presented in |Glampedakis et al.| 
|20l4| . 


3.7 Boundaries and boundary conditions 


For specifying the various regions of the neutron star, as well 
as the boundary conditions, we use the same assumptions 
as in Lander (20141, modified in a way that entrainment is 


included. 


3.7.1 The crust-core and surface boundaries 

The surface of the neutron star is defined by the contour of 
vanishing proton density 

Pp urf (w,z) = 0, (33) 

with a corresponding radius r e := r% q in the equatorial 
plane. The core of the neutron star extends from the center 
out to the proton density contour that is at a radius 0.9 rj? q 
in the equatorial plane, which also defines the crust-core 
boundary 

Pp c (w,a) := p p (0.9rP q ,0). (34) 

The crust is the region extending between the above crust- 
core boundary and the surface of the star. 


in such a way that the polynomial values coincide with the 
numerical ones at the pole and equator. Since we chose a 
second degree polynomial, a third point is needed as well, 
which we choose to be at the middle of the 9 direction (9 = 
7t/ 4). Therefore, we have 

co = C’ e , (38a) 


Cl 


Bcc - CO 
eq j 

a i M 


(38b) 


C2 = 


S mid „ „ „.mid 

cc C 0 C\U cc 

u£ id (w£ id - «c?) 


(38c) 


This polynomial approximation produces acceptable results, 
since it only induces a negligible error on the crust-core 
boundary. Since D p is a function of proton and ne utro n den¬ 
sities (D p = D p (pp, p n )) only, = 0 and hence (361 takes 
the following form 


y(u) = Dp C B cc (u) + 47t 


^crust 

p p 


Mn(«), 


(39) 


which relates y(u) with Mn(u). 

The second boundary condition ensures that B$ is con¬ 


tinuous at the crust-core boundary. Using (251 and (311 we 
obtain 


/ac(w) — (Hcl] c 


/n(u) 

B cc (u) ’ 


(40) 


which relates the superconducting and normal matter / 
functions. It is obvious that the independent functions are 


now two, instead of four. Even though (221 does not depend 
explicitly on the superconducting function M sc , this func¬ 
tion is used to obtain the equilibrium configurations, as will 
be shown later. The expression for M sc is found by solving 
(241 


M sc = i- (y(u) - BD P ], 


(41) 


and then substituting y from (391. 


3.7.2 The boundary conditions 


4 NUMERICAL METHOD 


The first boundary condition to be met is the continuity of 
the magnetic force on the crust-core boundary (we denote 
the crust-core boundary using the cc subscript) 


[,rVM BC ] cc =r t VM N ] cc . (35) 

Substituting (|24|) in (|35|) and assuming that Mn is a function 


of u only, we obtain (see A31 

[V (BDp )] cc = 


d Mi 


i ^crust 

— — 4n— -—— 

> core d u 


du 


Pp 


Vu 


(36) 


It is obvious that B cc = B cc (u). Since we do not know B as 
a function of u on the core-crust boundary explicitly, we will 
use a polynomial approximation for B cc , denoted by B cc ( u ), 
employing the exact same scheme as in Lander (2014]) 


B cc (u) = C 0 + ClU + C 2 u(u - Uec), 


(37) 


We use the Hachisu self-consistent field (HSCF) method 
(Hachisu|1986 I as presented by Lander et al. ( 2012|);|Lander 


(2014) with new convergence criteria to specify the equi¬ 
librium. First, we demonstrate these criteria and compare 
them to the original ones used by Hachisu (1986). Then, we 
describe the non-dimensional units, the plan of the method 
and finally the numerical implementation. The numerical 
code used to obtain the following solutions is an extension 
of the code presented in Palapanidis (2014). 


4.1 New convergence criteria 

In the style of the original HSCF method, we would require 
that in order for a numerical solution to be considered as a 
converged solution, the quantities 


where Co, Ci, c 2 are constants and is the equatorial value 
of u on the crust-core boundary. The constants are chosen 


A/Ip 


Pp max,new p-p max,old 

max,new 


(42) 
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A/I n — 


ACp = 


AC d if = 


Pn max,new Pn max,old 


P 


n max,new 


“P 


c D 


-c, 


p,old 


c D 


Cdif, 


■ — Cdif,( 


Cdif ,11 


(43) 


(44) 


(45) 


should all have become less than a specified value (usually 
of order 10 -6 ). Here, subscripts new and old denote consec¬ 
utive iterative steps. However, these criteria proved to be 
adequate only for relatively simple models, which converge 
to high accuracy after a small number of iterations, such as 


(2005). 


those presented in Hachisu (19861 and Tomimura & Eriguchi 


In the case of entrainment, the system of equations is 
more complex, requiring a larger number of iterations to 
achieve convergence. In order to ensure that all quantities 
relevant for the solution have converged, we introduce three 
new convergence quantities related to the proton density p p , 
the neutron density p n and the magnetic function u. The 
above physical quantities are fundamental, in the sense that 
all other quantities are functions of them. Motivated by the 
definition of the usual standard deviation, we introduce the 
following normalized error measures 


'Pp 


\ 


KDIV NDIV 


Y Y (* i,j new Pp i,j old ) 

i=l j=1 

NDIV ■ KDIV 


\ 


KDIV NDIV 


Y Y (Pn y new Pn i,j old ) 

*=1 3 = 1 _ 

NDIV • KDIV 


\ 


KDIV NDIV 


^ ^ ( (tq.j new 'U'ijj old) 

i=l j= 1 

NDIV ■ KDIV 


(46) 


(47) 


(48) 


where NDIV and KDIV are the number of equidistant grid 
points in the r and p := cos 9 directions respectively, while 
subscripts i, j specify the grid point. All of the above quan¬ 
tities should converge to zero. The first two are sensitive to 
the convergence of the fluid properties, while the last one is 
sensitive to the convergence of the magnetic field. 

When the magnetic field is not too strong, the density- 
dependent quantities (|46|) and (|47|) converge much faster 


than (481 which is related to the magnetic field. In such 


cases, we decouple the matter fields from the magnetic field 
when ( |46[ ) and ( |47| ) have either reaches a plateau or have 
become sufficiently small and continue iterating only the 
magnetic field equation, until (481 also reach a plateau or 


becomes sufficiently small. The decoupling also prevents the 
growth of inherent numerical inaccuracies due to the finite 


precision of the method. We note that (461,(471 and (481 are 


global error indicators and are effectively averaged over the 
whole numerical grid. In contrast, when using the original 


HSCF-style criteria (421-(451, two of these quantities are 
local quantities and do not show a monotonous variation 
during convergence - see Section|5.2| Using the new criteria 


(461,(471 and (481 it is more straightforward to decouple the 


iteration of the magnetic field (provided the magnetic field 
is weak) from the iteration of the fluid quantities, once the 
latter have become of sufficient accuracy. 

4.2 Non-dimensional Units 

In order to solve the equations numerically, we choose a sys¬ 
tem of non-dimensional units, which is based on the max¬ 
imum density p max , the gravitational constant G and the 
(proton) equatorial radius r fq. In this system, the length, 
mass and time units are 

[L] = r e p q , 


[At] = (r p q ) 3 p n 


[T} = 


(49) 

(50) 


(51) 


\J~VkG(> max 

In Appendix |B2| we derive the dimensionless form of various 
quantities which are shown with “ “ ”. Using the dimension¬ 
less units the form of the equations does not change. Divid¬ 
ing the proton chemical potential with the maximum value 


/j max and using (15a), (15b) as well as the definitions for the 


dimensionless densities (Bill and (B12 1 , one obtains 


Pp — *p(0) 


Pp 


Up i 


(52) 


and similarly 


Pn = (1 - Zp(0)) pN- • (53) 

where £ p (0) is the proton fraction (* p := p p /p) at the center 
of the star. 


4.3 Plan of the method 

We specify the polytropic indices N p , N n and the ratio of 
polar to equatorial radii for protons f pol /feq, which is simply 
equal to f p ol , since the dimensionless value of ffq = 1. The 
equatorial radius of neutrons r" q is determined through the 


definition of the crust-core surface (34) and its dimensionless 
form is the ratio fe q /r p q . The functions Mn(m) and /n(u) 
and the superconductivity parameter h c , as well as the cen¬ 
tral proton fraction * p (0) and entrainment constant <5* are 
specified. For numerical reasons, we usually need to spec¬ 
ify an under-relaxation parameter u> < 1, to achieve conver¬ 
gence. The proton chemical potential vanishes at the surface 
of the star, while the neutron chem ical pote ntial vanishes at 
the crust-core surface. Using (101 and (111 we evaluate p p 
at f p q , f p ol (where £ P (f p q ,0) = /x p (1,0): = 0, p P (f p ol , 1) = 0 
respectively) and p n at fg q (where /t n (feq,0) = 0) and we 
derive the equations for the integration constants C p and 
Cdif 


C p = 4> g (l, 0) — M(l, 0), 

Cdif = Pp (f*eq , 0) - M(f“ q , 1) + N(N q , 1). 
The main iteration algorithm is: 

Main iteration 


(54) 

(55) 
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(i) Assign initial values 6 P = 1, p n = 1 and u = 1. 

(ii) Compute 4> g from (mil. 


(iii) Calculate II from u (271. 


(iv) Compute intermediate u from (B31 (i.e. Hint). Before 
evaluating, divide by II max and multiply again after integra¬ 
tion. 

(v) Employ under-relaxation to find the new u 


= (1 — cu)uint + CUUold, 


(56) 


(vi) Evaluate the proton integration constant C p from 
(541. 


(vii) Use the proton integral equation (101 to evaluate jl p . 
(viii) Evaluate the difference integration constant Can 


from (551. 


(ix) Using the difference integral equation (111, compute 


Pn 


(x) Compute new proton and neutron densities from (521 
and (531. 


(xi) Return to first step and use for the next iteration the 
new values of p p , p n and u, until er Pp , <r Pn converge. 

(xii) Calculate II from u (271. (In the case of decoupling, 


hereafter we use fixed p p and p n values.) 


(xiii) Compute intermediate u from (B31. Before evaluat¬ 
ing, divide by n max and multiply again after integration. 

(xiv) Employ under-relaxation to find the new u. 

(xv) Return to step (xii) until a„ has converged. 

We use the aforementioned scheme to solve the system of 
equations numerically, employing a two-dimensional grid in 
spherical polar coordinates, with points equidistant in r and 
p := cos 9. Solving the integral equation for the gravitational 
potential is already well known ( jHachisui : 1986[ ). The equa¬ 
tion for obtaining u as well as its numerical implementation 
are given in Appendix |B1| 


4.4 Virial test 


The scalar virial theorem derivation for the two-fluid model 
is similar to the single fluid (see[Chandrasekhar 196l) [Collins| 
19781 and is given by 


1 A 2 1 
2iP 


= 2 T + £, 


^ g + W + 3^ + 3^+£ Fn , (57) 


where I is the moment of inertia, T is the kinetic energy, 
£ma g is the magnetic energy given by 


^■mag 


r ■ F mag dV) 


(58) 


W is the gravitational energy of the system, U n and U p is 
the internal energy of neutrons and protons respectively and 
finally £f„ is the energy due to the presence of entrainment 
given by 

£f„ = f r • F n dV, (59) 

J star 


which, hereafter, we will refer to as the ’’entrainment en¬ 
ergy” . The form of the last equation is analogous to (581 and 


given by the derivation of the scalar virial theorem (Collins 
p78l ). 


Since the moment of inertia needs to be constant with 
respect to time for stationary solutions, the right side of the 
virial theorem vanishes. Notice that we consider only non¬ 
rotating models. Numerically, the virial theorem will differ 


from zero and so we construct the virial test quantity, 

£mag + W + 3+ 3 ^ + £ Fji 

W\ 

which provides a test for the global convergence of the algo¬ 
rithm. Notice that for magnetized stars, this is a meaningful 
test only when the magnetic energy is a measurable fraction 
of the total energy. 


VC = 


( 60 ) 


5 RESULTS 
5.1 Assumptions 

In order to obtain a numerical solution of an equilibrium 
model, we have to specify the functions and the various pa¬ 
rameters governing the system of equations. For the normal 
matter magnetic functions /n ( u ) and Mn (m) we have chosen 

/n(w) = a (u - Mint) C + 1 B (u - Mint) , (61) 

Mn(u) = ku, (62) 

where a, £ and k are constants, Mint corresponds to the value 
of the largest constant-M line that closes inside the star (i.e. 
the u value on the equator) and 0(u) is the Heaviside step 
function. We consider purely poloidal held models. The pa¬ 
rameter a determines the strength of the toroidal part of the 
magnetic held while the value of k is related to the strength 
of both the poloidal and toroidal part of the magnetic held 
in the mixed held case. In the specihc models shown below, 
we omitted the toroidal component, since its contribution is 


of the neutron star is chosen as M = 1.64 Mq for all config¬ 
urations studied here. 

Our numerical results are converted from non- 
dimensional units to Gaussian-CGS units using p max = 
10 15 gcm~ 3 and r)’ q = 15 km. The value of the gravita¬ 
tional constant is G = 6.673 x 10 _s cm 3 g _1 s -2 . The nu¬ 
merical grid consisted of 480 x 480 points, incorporating one 
quadrant (using equatorial symmetry and axisymmetry) and 
an underrelaxation parameter uo = 0.02 was necessary for 
achieving convergence. The initial guesses for all the grid 
points are p p = 1, p n = 1 and u = 1. 


typically very small compared to the poloidal held, see (Lan¬ 


der & Jones 2009 Lander 20141. We choose aj p (0) = 0.15, 


h c = 0.1 and polytropic indices N n = 1, N p = 1. The mass 


5.2 The poloidal held configuration with 
entrainment 

We focus on three different models with magnetic held 
strength of order 10 lo G at the pole, constructed with en¬ 
trainment constant of J* = 0.8,0.9 and 1.0 (the latter value 
corresponds to the entrainment-free case). Since 5* is related 
to the proton density, we chose these representing values to 
obtain the equilibria ( |Chamel] [2008]). Table [T| summarizes 
the main properties of the models. The energy related to 
the entrainment is only a small fraction compared to the 
magnetic energy. The magnetic energy for the entrainment 
and entrainment-free cases is comparable for equilibria with 
similar magnetic helds. 

Fig0 shows the magnetic held lines (actually, contours 
of the magnetic function u) for the above three models. The 
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Table 1. Properties of three models with entrainment (first two lines) and without (third line). The values of the convergence quantities 
log cr Pn , log c p p are shown at the last fully-coupled iteration while log cr u is shown at the end of all iterations. 


Model 

Si, 

£mag/ 1 W\ 

£*J\W\ 

|W| 

Spole (10 15 G) 

log cr u 

!og CTp p 

lo S 

* 

Virial test 

with entrainment 

0.8 

3.33E-05 

4.66E-07 

6.12E-02 

1.16 

-6.93 

-6.65 

-6.18 

0.0320 

1.48E-05 

with entrainment 

0.9 

1.37E-05 

3.42E-07 

6.11E-02 

1.21 

-6.88 

-5.96 

-5.71 

0.0280 

5.32E-06 

no entrainment 

1.0 

1.27E-05 

- 

6.11E-02 

1.34 

-6.58 

-7.75 

-8.27 

0.0278 

3.53E-06 
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Figure 1 . Poloidal magnetic field lines (contours of u ) for the magnetized two-fluid type-II superconducting neutron stars of Table |T| 
The left-most and middle model have entrainment constant of <5* = 0.8,0.9 respectively. The right-most is without entrainment. The 
black circle represents the surface of the star while the white circle represents the crust-core boundary. The color scale shows magnitude 
of the magnetic field in dimensionless units. The three cases differ mainly by the presence of kinks in the contours of u at the crust-core 
boundary. 



Figure 2. Comparison of magnetic field lines (contours of u) for 
the models with entrainment parameter 5* = 0.8 (black lines) 
and entrainment-free model (red lines). The blue line represents 
the crust-core boundary. 


The low value of the magnetic energy with respect to 
fluid energies in the above models, implies that after the fluid 
quantities have numerically converged to their most accurate 
values, one could consider them as fixed and continue the it¬ 
erations for the magnetic field only, in order to improve its 
numerical convergence. The convergence test quantities as 
a function of the iteration number are shown in Fig. [3] (left 
panel) for the model with 6 * = 0.8. After nearly 500 itera¬ 
tions, the accuracy of the fluid quantities does not improve 
anymore and we continue iterating only the magnetic field 
equation. It takes about 200 additional iterations for the 
magnetic field test quantity (blue line) to reach a plateau. In 
contrast, if one would continue to update the fluid quantities 
in the iteration process, this would only accumulate numer¬ 
ical truncation errors, which would not allow the magnetic 
field to converge in a satisfactory manner. 

For comparison, we provide a similar graph of the 
Hachisu-style convergence test quantities in Fig. E (right 
panel). It is obvious that these quantities vary significantly 
during iterations, which might result in erroneously termi¬ 
nating the iterations earlier than desired. 


main qualitative difference among these models is the pres¬ 
ence of kinks in the contours of u at the crust-core bound¬ 
ary, which are especially evident in the 5* = 0.8 model. 
These kinks emerge because of the finite precision of the 
numerical treatment of the crust-core boundary. Matching 
the right-hand side values of [22] and [30] on the crust-core 
boundary is less accurate than in the entrainment-free case. 
This happens because the type-II superconducting Grad- 
Shafranov equation (221 obtains a much more complicated 
form (due to non-constant entrainment parameter e*) than 
that in the entrainment-free case presented in Lander (20141. 
Fig ,[2]compares directly the w-contours between the 5* = 0.8 
model and the entrainment-free case. The largest differences 
occur near the crust-core boundary and in the crust region. 


5.3 Stronger poloidal field results 

Here we discuss models with a stronger magnetic field (3 to 
4 times larger at the pole than the models in the previous 
Section). We focus on two particular models, one with en¬ 
trainment 5^ = 0.8 and one without entrainment. The main 
properties for both equilibria are shown in Table [2] In both 
models, the magnetic energy and the entrainment energy are 
higher when compared to the results of Table [T] 

In both cases, the maximum value of the magnetic field 
in the converged solution does not appear at the center of 
the star, but is displaced along the z-axis to a point between 
0.3 R E q and 0.4 R eq (Fig.pJ). This off-center maximum occurs 
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Figure 3. Left panel: The new convergence test quantities, a u (blue), a Pjl (magenta) and a Pp (green) as a function of iteration number 
for the case <5* = 0.8. The density-related test quantities reach a plateau before 500 iterations. At that point, we decouple the iteration 
of the magnetic field, in order to improve its accuracy. It takes about 200 additional iterations for the magnetic field test quantity (blue 
line) to reach a plateau. Right panel: The Hachisu-style convergence test quantities, A/7 n (blue), ACdif (magenta), AC P (green), and 
A/Ip (red). Their non-monotonous behaviour can terminate the iteration process too early. 



X 


Figure 4. Poloidal field for a magnetized two fluid type-II super¬ 
conducting neutron star with entrainment (<5* = 0.8). The black 
lines represent the u contours. The black circular line represents 
the surface of the star while the white line represents the crust- 
core boundary. We show the magnitude of the magnetic field with 
color. The magnetic field on the pole is 3.73 X 10 15 G. Here, the 
u kinks on the crust core boundary are much more obvious than 
in the weaker field case. 


when .Bpoie exceed roughly 3 x 10 15 G. In addition, in these 
models with stronger magnetic field, the kinks in the u con¬ 
tours on the crust-core boundary are much more apparent. 

Another aspect of the stronger field is that the conver¬ 
gence error indicators obtain a plateau at somewhat larger 
values than for the models in the previous Section. This im¬ 
plies that for models with stronger magnetic fields one would 
need better resolution to achieve the same accuracy. 


6 DISCUSSION 

We constructed Newtonian equilibrium models for magne¬ 
tized, axisymetric neutron stars taking, into account type-II 


superconductivity of protons, superfluidity of neutrons and 
entrainment. We employed the MHD equations presented in 


Glampedakis et al. (20111, a set of equations accounting for 


the forces related to fluxtubes and vortices. As of today, this 
is the most “complete” set of type-II superconducting MHD 
equations in the Newtonian framework. The models pre¬ 
sented here are derived from these equations by neglecting 
the effect of rotation on the structure of the star. Thus, the 
numerical models presented here are “complete” in the sense 
that they employ all terms discussed in Glampedakis et al. 
(20111 for non-rotating fluids in equilibrium. As pointed out 
in the introduction this paper provides an additional step 
towards our understanding of neutron star equilibria. The 
physics of entrainment is likely be of importance in the case 
of fast rotating stars. 


The poloidal field is qualitatively similar to the 
entrainment-free cases. The energy due to entrainment is 
only a very small fraction of the total energy - about an 
order of magnitude less than the magnetic energy. The new 
convergence test criteria we introduce, allow for a more ac¬ 
curate magnetic field configuration to be obtained. The main 
difference between models with and without entrainment 
and magnetic field at the pole less than ~ 3 x 10 15 G is struc¬ 
ture of magnetic field lines along the crust-core boundary. In 
the entrainment case, kinks are present that do not appear 
in the entrainment-free case. Though, as previously men- 
tionted, these kinks are of numerical nature. For stronger 
field magnetic fields, these kinks are stronger while the mag¬ 
netic field maximum appears off-center, along the symmetry 
axis. 


It is possible to include rotation in our models, although 
even the simplest case, rigid rotation, requires a much more 
complicated form of equations § and Q as various terms 
related to vorticity emerge. This constitutes a future step 
for this model to be achieved as well as including realistic 
equations of state rather than polytropes. Then, instead of 
approximating <5* as a constant, a relation between <5* and 
p, as shown in Charnel (20081, can be employed as well. Fur- 
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Table 2. Results for the stronger purely poloidal field configuration with entrainment and without. The values of the convergence 
quantities log a Pp , log a Pp are shown at the last fully-coupled iteration while log a u is shown at the end of all iterations. 


Model 

6* 

fmag/ |W| 

£fJ\W\ 

|W't 

M 

Bpola (10 15 G) 

log <J U 

!°g “"Pp 

log £7p n 

K 

Virial test 

with entrainment 

0.8 

9.16E-05 

4.65E-06 

6.11E-02 

0.968 

3.73 

-5.90 

-5.96 

-5.63 

0.0490 

3.06E-05 

no entrainment 

1.0 

4.53E-05 

- 

6.12E-02 

0.968 

5.16 

-6.11 

-7.14 

-7.75 

0.0420 

4.31E-06 


thermore, a more detailed description of the different physi¬ 
cal properties of the inner and outer parts of the crust could 
be used as well. In the case of strong crustal entrainment as 
discussed in Charnel (20121 different configurations might 
arise. Finally, it is possible to model a magnetosphere by 
extending the magnetic function /w in the exterior of the 
neutron star, as shown in Glampedakis et al. (20141. 
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APPENDIX A: MATHEMATICAL 
DERIVATIONS 


(20a I to find that: 


1 

47T 


+B x (Wei xb)+ p p V (BD P ) ]. 


(Al) 


Following 


Lander 


j := V X B = 


12014), we define a unit current j as 
x + j^e^. 




(A2) 


Using ( |A2[ ) in (JaTJ and using vector calculus identities we 
obtain 


— 4-7T F, 


mag = PpV {BDp) - —J (V (zuB^Hd'j X Vu) 


+— 

w 




+ — 


1 ( VRd ■ Vr 


w \ vjB 


- H cl j$ Vu. (A3) 


Since for the magnetic force, F mag = p p VM sc (see Eq.j8| the 
right side of the aforementioned equation should not contain 
any components in the (/-direction. Thus, the toroidal part 


of (A31 should be zero 
V (zuB^Hd'j x Vu = 0. 


(A4) 


Setting now / 3C := zuB^Hd we obtain that / sc is a function 
of u (i.e. /sc = /sc(w)). We also set y := A-kM sc + BD P and 
from (|A3[) we obtain 


- Vy 


ppzu au zo \ ZDrSpp 


^)vu. (AS) 


It is evident if we take the cross product of both parts of ( |A5| ) 
with Vu that y = y(u). Substituting B$ from the definition 
of /(u) we obtain 


dy 

du ppTZ7// c l dw 


Bj,f s c d / 3[ 


ygcl-Va I Hal ~A ± 
rr-2 Pi n_ m z>.. J 4* ‘ 


ZU 2 Bpp 


^Pp ° 


(A6) 


Using now the definition (A2 I of the unit current j and 


(131 we obtain the relation between j pi and u 


n = - 


zuB 


A *u — — VS • Vu 
B 


(A7) 


Substituting the aforementioned equation in (A6), we ob¬ 
tain the type-II superconducting Grad-Shafranov equation, 
similar to the one presented in Lander (2014): 


Al Derivation of the Grad-Shafranov equation 

In order to derive the Grad-Shafranov equation for a type-II 
superconductor we substitute H c i = H c iB in (171 and use 


A _Vn-Vu _ 2 ^dy n2( d/sc 

A *U- n w Pp n du n /-d„> 


where n is defined through (271. 


(AS) 
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A2 Entrainment quantities 

The explicit relations for D p , D n as functions of p p and p n 
are 

n _ dH cl _ ; p p (5* — l) 2 + 2p p p n (S* — 1) <5* + p 2 <5* f A nA 
p d Pp c ' - j 


n - _ / 

-On •— — h c 

dp n 


{pp ($* ~ 1) + PnS*) 

Pi (6+ - l) 2 

Op 0* - 1) + PnS*) 2 


(A10) 


It is obvious that when <5* = 1 (i.e. e* = 1) D p = 1 and 
D n = 0. 

A3 The boundary conditions 


The first boundary condition (361 dictates that [BD p ] cc is 


proton densities only, it holds that 


d(BD p ) 


drt 


D ^ 
p d u 


(All) 


APPENDIX B: NUMERICS 

B1 The numerical implementation of the 
Grad-Shafranov equation 

The Grad-Shafranov equation can be written as a Poisson 


equation using (231 and hence it can be inverted and solved 
for u as 


u{r, p) = 


l 


F( r') 


47r sin <j> J all space | r' * 
where F( r) is dehned through 

( -T+^Prnf + £/&, 

F =\ 4 + 

[ 0, 


sin <j>'dV', 


core 

crust 

exterior 


(Bl) 


(B2) 


(note that the minus sign by inverting Poisson’s equation 
has been absorbed in F). The primed quantities describe 
source points while non-primed points are where u is cal¬ 
culated. Expanding in (r, p) coordinates and due to 


axisymmetry, (B31 takes the following form: 

+oo 

^2 fin-i(r',r) 


r+oo /*1 

u(r, At) = vj J dr' / dp 
J r' =0 J pi' =0 


l 


2 n (2 n — 1) 


P^n-l (A 4 ) P^n-l (A 4 ) 


F(r',p'), 


(B3) 


where P{ n {p) are the associated Legendre polynomials and 
f n {r',r) is the radial part of the expansion, given by 


fn(r',r) = 


r /n -\-1 5 


r > r 
r <r' 


(B4) 


We now discretize (B31 in order to obtain a relation 


that can be implemented numerically. Following the HSCF 


method (Hachisu 1986) we employ a 2D r vs. p grid with 


NDIVxKDIV points in the respective directions. The coor¬ 
dinate p ranges between 0 and 1 while the radial coordinate 
r ranges between 0 and r max (we set r max =16/15 to cover 


the whole star and a small exterior region). The grid points 
are thus given by 

U ~ 1) 


U r max NDTV _ i 

and 


Vi — 


(i-i) 


(B5) 


(B6) 


KDIV - 1 

We compute n = LMAX terms of the Legendre and associ¬ 
ated Legendre polynomials (here we choose LMAX = 16). 
Using the same notation as before, integrating (Bit using 
Simpon’s rule, over p' yields 


U, 


(i) 


KDIV-2 

E 


i 


3(KDIV - 1) 


{P2n-l(Pl)Fk,l 


+4 P2n — l{Pl+l)Fk,l+l + P2n-l{Pl+2)Fk,l+2) 


(BT) 


while integration over r' gives 
NDIV -2 

ur 


( 2 ) _ 

n ’ j ^ 3 {NDIV - 1) 


= E 


(rlf2n.-i{r k ,r j )Ul 1 J l 


+^rk +1 f2n-i{r k+1 ,r j )U^ ln 
+rl+2f2u-i(r k+ 2,r j )ul 1 l 2 . 
Then, u is given by 


(B8) 


Uij = r. 


_ LM A.Ji 

E 


1 -r(2) 


K, P 2n -l(lM), (B9) 


^ 2n(2n-l) 

R=1 

where U)p n and U^ 2 j are intermediate integration quantities. 

B2 Dimensionless quantities 

Here we list various quantities in their dimensionless form 
r=^p~, (BIO) 

Pp = i (Bll) 

(B12) 
(B13) 
(B14) 
(B15) 
(B16) 
(B17) 
(B18) 


Pp 

P max 
Pn 


<i> = __ 

47tG (rf q ) 2 p max ’ 

A C p 


Cdif = 


47tG (rfq) p n 

Gdif 


p \2 ’ 

“ i n 

max 


Pp 


Pn — 


47tG (rfq) p , 

_ Pp _ 

47tG (rl'q) 2 p max ’ 

Pn 

4ttG (rl q) 2 p max ’ 
u 

^eq) Pmax 
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M = 



K 



a = 


a 

1 ’ 


(B19) 

(B20) 


(B21) 


y/ 47rGp m ax ( r eq) 

while for all quantities with energy units the dimensionless 
form is 

E 


E = 


47 tG (rl’q) pn 


(B22) 


This paper has been typeset from a T1]X/ LMJrjX hie prepared 
by the author. 



